!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine roe_x (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                           &
   ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, ny, gm)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny
   real(8), intent(in) :: gm
   real(8), dimension(nx, ny), intent(in) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l
   real(8), dimension(nx, ny), intent(in) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r
   real(8), dimension(nx, ny), intent(inout) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen

   integer(4) :: i, ip1, j, i_dim
   real(8), dimension(nx, ny) :: ro_a, mx_a, my_a, mz_a, bx_a, by_a, bz_a, en_a
   real(8), dimension(8) :: var_c
   real(8), dimension(8, nx, ny) :: eigenvalue
   real(8), dimension(nx, ny) :: fro_l, fmx_l, fmy_l, fmz_l, fbx_l, fby_l, fbz_l, fen_l
   real(8), dimension(nx, ny) :: fro_r, fmx_r, fmy_r, fmz_r, fbx_r, fby_r, fbz_r, fen_r
   real(8), dimension(8, nx, ny) :: var_d
   real(8), dimension(8, 8, nx, ny) :: eigenvector_l, eigenvector_r
   real(8) :: gmm1, gmm2
   real(8) :: vx_l, vy_l, vz_l, b2_l, v2_l, pr_l, h_l, ro_sl
   real(8) :: vx_r, vy_r, vz_r, b2_r, v2_r, pr_r, h_r, ro_sr

!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0
   gmm2 = gm - 2.0d0

   do j = 1, ny
   do i = 1, nx - 1

      ip1 = i + 1

      vx_l = mx_l(i, j) / ro_l(i, j)
      vy_l = my_l(i, j) / ro_l(i, j)
      vz_l = mz_l(i, j) / ro_l(i, j)
      b2_l = bx_l(i, j) * bx_l(i, j) + by_l(i, j) * by_l(i, j) + bz_l(i, j) * bz_l(i, j)
      v2_l = vx_l * vx_l + vy_l * vy_l + vz_l * vz_l
      pr_l = (en_l(i, j) - v2_l * ro_l(i, j) / 2.0d0 - b2_l / 2.0d0) * gmm1
      h_l = (ro_l(i, j) * v2_l / 2.0d0 + gm * pr_l / gmm1 + b2_l) / ro_l(i, j)
      fro_l(i, j) = mx_l(i, j)
	   fmx_l(i, j) = vx_l * mx_l(i, j) + pr_l + b2_l / 2.0d0 - bx_l(i, j) * bx_l(i, j)
	   fmy_l(i, j) = vx_l * my_l(i, j) - bx_l(i, j) * by_l(i, j)
	   fmz_l(i, j) = vx_l * mz_l(i, j) - bx_l(i, j) * bz_l(i, j)
	   fbx_l(i, j) = 0.0d0
	   fby_l(i, j) = vx_l * by_l(i, j) - bx_l(i, j) * vy_l
	   fbz_l(i, j) = vx_l * bz_l(i, j) - bx_l(i, j) * vz_l
	   fen_l(i, j) = (en_l(i, j) + pr_l + b2_l / 2.0d0) * vx_l - (bx_l(i, j) * vx_l + by_l(i, j) * vy_l + bz_l(i, j) * vz_l) * bx_l(i, j)

      vx_r = mx_r(ip1, j) / ro_r(ip1, j)
      vy_r = my_r(ip1, j) / ro_r(ip1, j)
      vz_r = mz_r(ip1, j) / ro_r(ip1, j)
      b2_r = bx_r(ip1, j) * bx_r(ip1, j) + by_r(ip1, j) * by_r(ip1, j) + bz_r(ip1, j) * bz_r(ip1, j)
      v2_r = vx_r * vx_r + vy_r * vy_r + vz_r * vz_r
      pr_r = (en_r(ip1, j) - v2_r * ro_r(ip1, j) / 2.0d0 - b2_r / 2.0d0) * gmm1
      h_r = (ro_r(i, j) * v2_r / 2.0d0 + gm * pr_r / gmm1 + b2_r) / ro_r(ip1, j)
      fro_r(i, j) = mx_r(ip1, j)
	   fmx_r(i, j) = vx_r * mx_r(ip1, j) + pr_r + b2_r / 2.0d0 - bx_r(ip1, j) * bx_r(ip1, j)
	   fmy_r(i, j) = vx_r * my_r(ip1, j) - bx_r(ip1, j) * by_r(ip1, j)
	   fmz_r(i, j) = vx_r * mz_r(ip1, j) - bx_r(ip1, j) * bz_r(ip1, j)
	   fbx_r(i, j) = 0.0d0
	   fby_r(i, j) = vx_r * by_r(ip1, j) - bx_r(ip1, j) * vy_r
	   fbz_r(i, j) = vx_r * bz_r(ip1, j) - bx_r(ip1, j) * vz_r
	   fen_r(i, j) = (en_r(ip1, j) + pr_r + b2_r / 2.0d0) * vx_r - (bx_r(ip1, j) * vx_r + by_r(ip1, j) * vy_r + bz_r(ip1, j) * vz_r) * bx_r(ip1, j)

!      ro_a(i, j) = (ro_l(i, j) + ro_r(ip1, j)) / 2.0d0
!      mx_a(i, j) = (mx_l(i, j) + mx_r(ip1, j)) / 2.0d0
!      my_a(i, j) = (my_l(i, j) + my_r(ip1, j)) / 2.0d0
!      mz_a(i, j) = (mz_l(i, j) + mz_r(ip1, j)) / 2.0d0
!      bx_a(i, j) = (bx_l(i, j) + bx_r(ip1, j)) / 2.0d0
!      by_a(i, j) = (by_l(i, j) + by_r(ip1, j)) / 2.0d0
!      bz_a(i, j) = (bz_l(i, j) + bz_r(ip1, j)) / 2.0d0
!      en_a(i, j) = (en_l(i, j) + en_r(ip1, j)) / 2.0d0

      ro_sl = sqrt (ro_l(i, j))
      ro_sr = sqrt (ro_r(ip1, j))

      ro_a(i, j) = (ro_sr * ro_l(i, j)  + ro_sl * ro_r(ip1, j)) / (ro_sl + ro_sr)
      mx_a(i, j) = (ro_sl * vx_l  + ro_sr * vx_r) / (ro_sl + ro_sr) * ro_a(i, j)
      my_a(i, j) = (ro_sl * vy_l  + ro_sr * vy_r) / (ro_sl + ro_sr) * ro_a(i, j)
      mz_a(i, j) = (ro_sl * vz_l  + ro_sr * vz_r) / (ro_sl + ro_sr) * ro_a(i, j)
      bx_a(i, j) = (ro_sr * bx_l(i, j)  + ro_sl * bx_r(ip1, j)) / (ro_sl + ro_sr)
      by_a(i, j) = (ro_sr * by_l(i, j)  + ro_sl * by_r(ip1, j)) / (ro_sl + ro_sr)
      bz_a(i, j) = (ro_sr * bz_l(i, j)  + ro_sl * bz_r(ip1, j)) / (ro_sl + ro_sr)
      en_a(i, j) = ((ro_sl * h_l  + ro_sr * h_r) / (ro_sl + ro_sr) * ro_a(i, j) +                                                  &
         gmm1 * (mx_a(i, j) * mx_a(i, j) + my_a(i, j) * my_a(i, j) + mz_a(i, j) * mz_a(i, j)) / ro_a(i, j) / 2.0d0 +               &
         gmm2 * (bx_a(i, j) * bx_a(i, j) + by_a(i, j) * by_a(i, j) + bz_a(i, j) * bz_a(i, j)) / 2.0d0) / gm

      var_d(1, i, j) = ro_l(i, j) - ro_r(ip1, j)
      var_d(2, i, j) = mx_l(i, j) - mx_r(ip1, j)
      var_d(3, i, j) = my_l(i, j) - my_r(ip1, j)
      var_d(4, i, j) = mz_l(i, j) - mz_r(ip1, j)
      var_d(5, i, j) = bx_l(i, j) - bx_r(ip1, j)
      var_d(6, i, j) = by_l(i, j) - by_r(ip1, j)
      var_d(7, i, j) = bz_l(i, j) - bz_r(ip1, j)
      var_d(8, i, j) = en_l(i, j) - en_r(ip1, j)
   enddo
   enddo

   call eigenspace (eigenvalue, eigenvector_l, eigenvector_r, ro_a, mx_a, my_a, mz_a, bx_a, by_a, bz_a, en_a, nx, ny, gm)

   do j = 1, ny
   do i = 1, nx
      do i_dim = 1, 8
         var_c(i_dim) = sum (eigenvector_l(:, i_dim, i, j) * var_d(:, i, j))
      enddo

      var_c = var_c * abs (eigenvalue(:, i, j))

      fro(i, j) = (fro_l(i, j) + fro_r(i, j) + sum (var_c * eigenvector_r(:, 1, i, j))) * 0.5d0
      fmx(i, j) = (fmx_l(i, j) + fmx_r(i, j) + sum (var_c * eigenvector_r(:, 2, i, j))) * 0.5d0
      fmy(i, j) = (fmy_l(i, j) + fmy_r(i, j) + sum (var_c * eigenvector_r(:, 3, i, j))) * 0.5d0
      fmz(i, j) = (fmz_l(i, j) + fmz_r(i, j) + sum (var_c * eigenvector_r(:, 4, i, j))) * 0.5d0
      fbx(i, j) = (fbx_l(i, j) + fbx_r(i, j) + sum (var_c * eigenvector_r(:, 5, i, j))) * 0.5d0
      fby(i, j) = (fby_l(i, j) + fby_r(i, j) + sum (var_c * eigenvector_r(:, 6, i, j))) * 0.5d0
      fbz(i, j) = (fbz_l(i, j) + fbz_r(i, j) + sum (var_c * eigenvector_r(:, 7, i, j))) * 0.5d0
      fen(i, j) = (fen_l(i, j) + fen_r(i, j) + sum (var_c * eigenvector_r(:, 8, i, j))) * 0.5d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine roe_x
